MODIS NRT Global Flood Product

Introduction

The MODIS/Aqua+Terra Global Flood Product L3 Near Real Time (NRT) 250m Global Flood Product (MCDWD_L3_NRT) (beta) provides daily maps of flooding globally. The product is provided over 3 compositing periods (1-day, 2-day, and 3-day) to minimize the impact of clouds and more rigorously identify flood water (the best composite will depend on the cloudiness for a particular event).

NASA EARTHDATA

CRM SEARCH

packages_to_check <- c("stars", "httr", "jsonlite", "tmap")

# Check and install packages
for (package_name in packages_to_check) {
  if (!package_name %in% rownames(installed.packages())) {
    install.packages(package_name)
    cat(paste("Package", package_name, "installed.\n"))
  } else {
    cat(paste("Package", package_name, "is already installed.\n"))
  }
  library(package_name, character.only = TRUE)
}

Package stars is already installed.

Loading required package: abind

Loading required package: sf

Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE

Package httr is already installed.
Package jsonlite is already installed.
Package tmap is already installed.

Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
remotes::install_github('r-tmap/tmap')

#in case tmap does not install
#remotes::install_github('r-tmap/tmap')

Reading the Data

Check what days are available: LANCE NRT

Based on availability, edit the year_day variable YYYY-DD. Example: ‘2022-01’

#add the year and date you want to search for (YYYY-DD, 2022-01)
year_day <- '2023-336'

Determine tiles of interest: MODIS NRT Tile Map

Based on availability, edit the tile_code variable:

#add tile code from the map above (written as h00v00)
tile_code <- 'h05v05'

This is the NRT Flood F3 (MCDWD_L3_F3) API URL:

API_URL <- paste0('https://nrt3.modaps.eosdis.nasa.gov/api/v2/content/details?products=MCDWD_L3_F3_NRT&archiveSets=61&temporalRanges=')

We can combine the API URL above with the year_day provided and print the available datasets:

#pasting together URL and year_day
url <- paste0(API_URL, year_day)
print(url)

[1] "https://nrt3.modaps.eosdis.nasa.gov/api/v2/content/details?products=MCDWD_L3_F3_NRT&archiveSets=61&temporalRanges=2023-336"

Load Data

Access the NASA Earthdata with the GET function:

# Make the GET request
response <- httr::GET(url)

Check response status from the GET function:

response

Response [https://nrt3.modaps.eosdis.nasa.gov/api/v2/content/details?products=MCDWD_L3_F3_NRT&archiveSets=61&temporalRanges=2023-336]
  Date: 2023-12-08 07:32
  Status: 200
  Content-Type: application/json;charset=UTF-8
  Size: 112 kB

# Check the response status
# Check the response status
if (http_status(response)$category == "Success") {
  # Parse the response JSON
  data <- content(response, as = "text", encoding = "UTF-8")
  data_parsed <- jsonlite::fromJSON(data)
  #filter for the tile code
  content_items <- data_parsed$content[grepl(tile_code, data_parsed$content$name, ignore.case = TRUE), ]
  #content_items <- data_parsed[grepl("h09v06", data_parsed$name, ignore.case = TRUE), ]
} else {
  print("Request failed with status code", http_status(response)$status_code)
}

print(content_items)

    archiveSets      cksum    collections               dataDay
116          61 1275757885 modis-nrt-c6.1 2023-336 = 2023-12-02
                                                                                          downloadsLink
116 https://nrt3.modaps.eosdis.nasa.gov/api/v2/content/archives/MCDWD_L3_F3_NRT.A2023336.h05v05.061.tif
        fileId                           md5sum      mtime
116 2051979420 64a8ba49c9893388c1da7217df15ec77 1701612811
                                       name        products resourceType
116 MCDWD_L3_F3_NRT.A2023336.h05v05.061.tif MCDWD_L3_F3_NRT         File
                                                               self    size
116 /api/v2/content/details/MCDWD_L3_F3_NRT.A2023336.h05v05.061.tif 1080978

Select the URL from the ‘downloadsLink’ column in the content_items list:

download_link <- content_items$downloadsLink
print(download_link)

[1] "https://nrt3.modaps.eosdis.nasa.gov/api/v2/content/archives/MCDWD_L3_F3_NRT.A2023336.h05v05.061.tif"

Use the “read_stars()” function from the “stars” R Library to read the geoTiff raster.
The raster is assigned to the “x” variable:

raster_df <- stars::read_stars(download_link)

Warning: ignoring unrecognized unit: none

my_raster <- st_set_crs(raster_df, st_crs("EPSG:4326"))

Warning in `st_crs<-.dimensions`(`*tmp*`, value = value): replacing CRS does
not reproject data: use st_transform, or st_warp to warp to a new CRS

st_crs(my_raster)

Coordinate Reference System:
  User input: EPSG:4326 
  wkt:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]

Visualizing NRT Flood Data

Plot the raster to quickly view it:

plot(my_raster, axes = TRUE)

downsample set to 3

Create NRT Flood Classification

Refer to the MODIS NRT Global Flood Product User Guide for more information.

NRT Flood data has 5 classifications:

Code	Definition
0	No Water
1	Surface Water
2	Recurring flood1
3	Flood (unusual)
255	Insufficient data

Create a classified legend; however, the NRT Flood data is stored in decimal numbers (aka floating-point). Create class breaks dividing the data by these breaks, and corresponding colors and labels:

class_breaks <- c(-Inf, 0.1, 1.1, 2.1, 3.1)
colors <- c( "gray", "blue", "yellow","red")

labels = c("0: No Water", "1: Surface Water", "2: Recurring flood", "3: Flood (unusual)")

Add a title for the plot that includes the year, day of year, and tile code:

title = paste("NRT Flood", year_day, tile_code)

Generate a plot from the tmap library using the tm_shape() function. With style as “cat,” meaning categorical. T

tmap_mode("view")

tmap mode set to interactive viewing

## tmap mode set to plotting
tm_shape(raster_df, style="cat" )+
  tm_raster(palette = c(colors), 
  title = title, 
  breaks = class_breaks,
  labels = labels )+
  tm_basemap(server = "Esri.WorldImagery") +
  tm_layout(legend.outside = TRUE)

stars object downsampled to 1000 by 1000 cells. See tm_shape manual (argument raster.downsample)

Footnotes

1. Value 2 (Recurring flood) is not populated in the beta release.